#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Appendix analysis
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

rm(list = ls())

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Packages, functions, parameters ####
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

## Packages ----

# packages: Core and functional packages
library(tidyverse)
library(lubridate)
library(zoo)

# packages: Matching packages
library(PanelMatch)

# packages: Mapping
library(choroplethrMaps)


## Functions ----

# function `trendplot_my': Plot for logged detainers (month-year version)
trendplot_my <- function(byvar, treat.set, control.set, labels, ylab) {
  
  # get treated means
  t <- data.frame(y = tapply(treat.set[,byvar], treat.set$t, mean, na.rm = T))
  t$x <- rownames(t)
  t$t <- "t"
  
  # get control means 
  c <- data.frame(y = tapply(control.set[,byvar], control.set$t, mean, na.rm = T))
  c$x <- rownames(c)
  c$t <- "c"
  
  # put them together
  toplot <- rbind(t, c)
  
  # plot
  toplot$t <- factor(toplot$t, levels = c("c", "t"))
  toplot$x <- as.numeric(as.character(toplot$x))
  
  g <- ggplot(data = toplot, aes(x = x, y = y, colour = t, group = t)) + 
    geom_point() + 
    geom_smooth(data = toplot[toplot$x<0,], method = "loess", se = FALSE) + 
    geom_smooth(data = toplot[toplot$x>=0,], method = "loess", se = FALSE) + 
    geom_vline(xintercept = 0) + 
    theme_minimal() + 
    xlab("Months to 287(g) signing") + 
    ylab(ylab) + 
    theme(legend.position="bottom") + 
    labs(color = "") +
    scale_color_manual(labels = labels,
                       values = c("c" = "gray30", "t" = "red3"))
  print(g)
}

# function `panel_match': Custom matching function 
panel_match <- function(byvar, lag = 12, nmatches = 10, lead = c(0:36), df, match.vars = "", 
                        time.id, unit.id, treat) {
  
  # make empty list
  matchlist <- list()
  
  # write matching formula
  match.vars <- c("", match.vars)
  covs.formula <- as.formula(paste0("~ I(lag(", byvar, ", 1:12))", paste(match.vars, collapse =" + ")))
  
  # get all treated units
  treat.ids <- unique(df[df[,treat]==1, unit.id])
  
  for (i in 1:length(treat.ids)) {
    # eliminate potential controls that would be treated over the analysis window
    time.treat <- min(unique(df$date1[df[,unit.id]==treat.ids[i]]), na.rm = TRUE)
    period.treat <- unique(df[df$my==as.yearmon(time.treat),time.id])
    start.period <- period.treat - lag
    end.period <- period.treat + max(lead)
    sub.treat <- df[df[,unit.id]==treat.ids[i] & df$period %in% c(start.period:end.period),]
    lose.ids <- unique(df[df[,treat]==1 & df[,time.id] %in% c(start.period:end.period), unit.id])
    sub.control <- df[!df[,unit.id] %in% lose.ids & df$period %in% c(start.period:end.period),]
    sub <- rbind(sub.treat, sub.control)
    
    # run matching
    cat(paste0("running matching on unit ", i, "\n"))
    matches.maha <- PanelMatch(lag = lag, 
                               lead = lead, 
                               time.id = time.id,
                               unit.id = unit.id, 
                               treatment = treat,
                               outcome.var = byvar,
                               refinement.method = "mahalanobis", 
                               covs.formula = covs.formula,
                               size.match = nmatches,
                               qoi = "att", 
                               data = sub)
    
    # save to list
    matchlist[[i]] <- matches.maha$att
  }
  
  # name and output list 
  names(matchlist) <- treat.ids
  return(matchlist)
}

# function `construct.sets`: use matches to construct datasets, aggregated to month
construct.sets <- function(mm, data) {
  treat.set.matched <- control.set.matched <- NULL
  matches <- names(mm)
  for (m in 1:length(matches)) {
    match <- mm[[m]]
    control.units <- as.numeric(as.vector(match[[1]][attributes(match[[1]])$weights>0]))
    treat.unit <- summary(match)$overview$id
    treat.period <- summary(match)$overview$period
    t0 <- as.numeric(as.character(treat.period))
    mindate <- t0-12
    maxdate <- t0+12
    sub.treat <- data[data$period>=mindate & data$period<=maxdate & data$id==treat.unit,]
    sub.control <- data[data$period>=mindate & data$period<=maxdate & data$id %in% control.units,]
    sub.treat$t <- sub.treat$period-t0
    sub.control$t <- sub.control$period-t0
    treat.set.matched <- rbind(treat.set.matched, sub.treat)
    control.set.matched <- rbind(control.set.matched, sub.control)
  }
  return(list(treat.set = treat.set.matched, control.set = control.set.matched))
}

# function `compute.sets.year`: use matches to construct datasets, aggregated to year
construct.sets.year <- function(mm, data_month, data_year, years_before, years_after) {
  treat.set.matched <- control.set.matched <- NULL
  matches <- names(mm)
  summary <- data.frame(matrix(NA, ncol = 3, nrow = length(matches)))
  names(summary) <- c("year", "treat", "matches")
  for (m in 1:length(matches)) {
    match <- mm[[m]]
    control.units <- as.numeric(as.vector(match[[1]][attributes(match[[1]])$weights>0]))
    treat.unit <- summary(match)$overview$id
    treat.period <- summary(match)$overview$period
    # convert them to format for by-year dataset
    treat.sc <- unique(data_month$sc[data_month$id==treat.unit])
    control.sc <- unique(data_month$sc[data_month$id %in% control.units])
    treat.year <- unique(data_month$year[data_month$period==treat.period])
    
    # construct summary dataset
    summary$year[m] <- treat.year
    summary$treat[m] <- as.character(treat.sc)
    summary$matches[m] <- paste(sort(control.sc), collapse = ", ")
    
    # construct analysis dataset
    t0 <- as.numeric(treat.year)
    mindate <- t0 - years_before
    maxdate <- t0 + years_after
    sub.treat <- data_year[data_year$year>=mindate & data_year$year<=maxdate & data_year$sc==treat.sc,]
    sub.control <- data_year[data_year$year>=mindate & data_year$year<=maxdate & data_year$sc %in% control.sc,]
    sub.treat$t <- sub.treat$year-t0
    sub.control$t <- sub.control$year-t0
    treat.set.matched <- rbind(treat.set.matched, sub.treat)
    control.set.matched <- rbind(control.set.matched, sub.control)
  }
  return(list(treat.set = treat.set.matched, control.set = control.set.matched, summary = summary))
}

# function `compute.did`: Difference-in-differences estimation
compute.did <- function(mm.obj, byvar, data, sets, time.id, unit.id, treat, lag, lead) {
  
  # manual diff-in-diff estimation
  did_est_t0 <- mean(sets$treat.set[sets$treat.set$t==0, byvar], na.rm = TRUE) - mean(sets$treat.set[sets$treat.set$t==-1, byvar], na.rm = TRUE) - 
    (mean(sets$control.set[sets$control.set$t==0, byvar], na.rm = TRUE) - mean(sets$control.set[sets$control.set$t==-1, byvar], na.rm = TRUE))
  
  did_est_t1 <- mean(sets$treat.set[sets$treat.set$t==1, byvar], na.rm = TRUE) - mean(sets$treat.set[sets$treat.set$t==-1, byvar], na.rm = TRUE) - 
    (mean(sets$control.set[sets$control.set$t==1, byvar], na.rm = TRUE) - mean(sets$control.set[sets$control.set$t==-1, byvar], na.rm = TRUE))
  
  did_est_t2 <- mean(sets$treat.set[sets$treat.set$t==2, byvar], na.rm = TRUE) - mean(sets$treat.set[sets$treat.set$t==-1, byvar], na.rm = TRUE) - 
    (mean(sets$control.set[sets$control.set$t==2, byvar], na.rm = TRUE) - mean(sets$control.set[sets$control.set$t==-1, byvar], na.rm = TRUE))
  
  did_est_t3 <- mean(sets$treat.set[sets$treat.set$t==3, byvar], na.rm = TRUE) - mean(sets$treat.set[sets$treat.set$t==-1, byvar], na.rm = TRUE) - 
    (mean(sets$control.set[sets$control.set$t==3, byvar], na.rm = TRUE) - mean(sets$control.set[sets$control.set$t==-1, byvar], na.rm = TRUE))
  
  manual_ests <- c(did_est_t0, did_est_t1, did_est_t2, did_est_t3)
  names(manual_ests) <- c("t=0", "t=1", "t=2", "t=3")
  
  # repeat with panelmatch to get standard errors, check 
  covs.formula <- as.formula(paste0("~ I(lag(", byvar, ", 1))")) # doesn't matter, will swap
  
  # run PanelMatch
  mm_temp <- PanelMatch(lag = lag,  
                        lead = lead, 
                        time.id = time.id,
                        unit.id = unit.id, 
                        treatment = treat,
                        outcome.var = byvar,
                        refinement.method = "mahalanobis", 
                        covs.formula = covs.formula,
                        qoi = "att", 
                        data = data)
  
  # swap in matches from before
  treats <- as.character(names(mm_temp$att))
  treat.ids <- str_split(treats, pattern = "\\.", simplify = TRUE)[,1]
  ids.in.match.obj <- names(mm.obj)
  for (m in 1:length(treats)) {
    weights_to_replace <- attributes(mm_temp$att[[m]])$weights
    ind.to.use <- which(ids.in.match.obj==treat.ids[m])
    weights_from_before <- attributes(mm.obj[[ind.to.use]][[1]])$weights
    for (w in c(1:length(weights_to_replace))) {
      weights_to_replace[w] <- weights_from_before[names(weights_to_replace[w])]
    }
    # cat(paste0("number of obs found = ", sum(weights_to_replace>0, na.rm = TRUE), "\n"))
    attributes(mm_temp$att[[m]])$weights <- weights_to_replace
  }
  
  pe <- PanelEstimate(sets = mm_temp, data = data, confidence.level = c(.95, .99, .999))
  pe <- data.frame(summary(pe)$summary)
  
  # compute star for table
  pe$star <- ifelse(sign(pe[,"X0.05."])==sign(pe[,"X99.95."]), "***", 
                    ifelse(sign(pe[,"X0.5."])==sign(pe[,"X99.5."]), "**",
                           ifelse(sign(pe[,"X2.5."])==sign(pe[,"X97.5."]), "*", "")))
  
  return(list(pe = pe, manual = manual_ests))
}

# statistical significance "stars"
star <- function(t, c, var) {
  star <- ifelse(t.test(t[t$t==-1, var], c[c$t==-1, var])$p.value<.001, "***", 
                 ifelse(t.test(t[t$t==-1, var], c[c$t==-1, var])$p.value>=.001 & t.test(t[t$t==-1, var], c[c$t==-1, var])$p.value<.1, "**", 
                        ifelse(t.test(t[t$t==-1, var], c[c$t==-1, var])$p.value>=.01 & t.test(t[t$t==-1, var], c[c$t==-1, var])$p.value<.05, "*", "")))
  return(star)
}


## Global parameters ----

# set seed
set.seed(123)

# set time zone (to avoid lubridate errors)
Sys.setenv(TZ = "America/New_York")

# define plot colors 
red_mit = '#A31F34'
red_light = '#A9606C'
blue_mit = '#315485'
grey_light= '#C2C0BF'
grey_dark = '#8A8B8C'
black = '#353132'

# define map colors
bloodred <- rgb(0.502,0.000,0.000, alpha=.75)
lightred <- rgb(0.502,0.000,0.000, alpha=.25)
lightgray <- rgb(0.8,0.8,0.8)



#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Load datasets  #### 
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

## Month-year level data set ----

# load month-year data 
full_my <- read.csv("full_my.csv")

# convert my back to the right format
full_my$my <- as.yearmon(full_my$my)
full_my$date1 <- ymd(full_my$date1)

# logged controls
full_my$pop_log <- log(full_my$pop + 1) 
full_my$total_crimes_log <- log(full_my$total_crimes + 1)

# crime rate per thousand citizens 
full_my$crime_rate <- full_my$total_crimes/(full_my$pop/1000)


## Month-year level subset for matching ----

# subset 
subset <- full_my[!is.na(full_my$treat),]

# create county-level estimate of undocumented by multiplying county-level foreign-born by state-level proportion of foreign-born who are undocumented
subset$undoc_county <- subset$foreign * subset$prop_undoc_fb

# detainers per thousand undocumented
subset$ndet_perundoc <- subset$ndetainees/(subset$undoc_county + 1)*1000

# detainers per thousand foreign-born 
subset$ndet_perforeign <- subset$ndetainees/(subset$foreign + 1)*1000

# logged detaienrs
subset$log_detainees <- log(subset$ndetainees + 1)

# get dataset of dates of switch 
dates <- unique(subset[!is.na(subset$date), c("sc", "date")])
dates$date <- ymd(dates$date)

# keep only first date for each county
dates <- dates %>% group_by(sc) %>% summarise(min_date = min(date))
dates$switch_mo <- as.yearmon(dates$min_date)
dates$year <- year(dates$min_date)

touse <- dates$sc

subset$id <- as.numeric(as.factor(subset$sc))
subset <- subset[order(subset$sc, subset$my),]
subset$period <- as.integer(as.factor(subset$my))

# change treatment here so that it's not aggregated by year
subset$date1_my <- as.yearmon(subset$date1)
subset$treat[subset$my<subset$date1_my] <- 0 


## Year level dataset ----

# load data 
full <- read.csv("full_y.csv") %>%
  select(
    year, GEO_id, date1, sc, denied, treat, control, signer,
    pop, pct_dem, prot2006, sec_comm,
    foreign, delta5_foreign, prop_undoc_fb,
    ndetainees, ndetainees_county, 
    ndetainees_federal, ndetainees_ice, ndetainees_local, ndetainees_state,
    ncharged, ncharged_county, ncharged_high_county, ncharged_med_county, ncharged_low_county
  )


# create county-level estimate of undocumented by multiplying
#   county-level foreign-born by state-level proportion of 
#   foreign-born who are undocumented
full$undoc_county <- full$foreign * full$prop_undoc_fb

### outcomes 

# logged detainer requests
full$log_detainees <- log(full$ndetainees + 1)
full$log_detainees_county <- log(full$ndetainees_county + 1)
full$log_detainees_noncounty <- log(full$ndetainees_federal + full$ndetainees_ice + full$ndetainees_local + full$ndetainees_state + 1)

# logged charged
full$log_notcharged_county <- log(full$ndetainees_county - full$ncharged_county + 1)
full$log_ncharged_county <- log(full$ncharged_county + 1)
full$log_ncharged_high_county <- log(full$ncharged_high_county + 1)
full$log_ncharged_med_county <- log(full$ncharged_med_county + 1)
full$log_ncharged_low_county <- log(full$ncharged_low_county + 1)

# get rid of observations with missing treatment status 
full <- full[!is.na(full$treat),]

# merge on id variable from state-year dataset
ids <- unique(subset[,c("sc", "id")])
full <- merge(full, ids, by = "sc", all.x = TRUE)

# create period variable
full$period <- as.integer(full$year - 2001)

# sort
full <- full[order(full$sc, full$period),]



#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Construct matched data sets  #### 
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


## Use PanelMatch to create matched sets ----

match.vars <- c("pop", "foreign")
mm3 <- panel_match(byvar = "log_detainees", treat = "treat", time.id = "period", unit.id = "id", nmatches = 3, match.vars = match.vars, df = subset)

# save matches
matches <- list(mm3 = mm3)


## Construct datasets based on matches ----

# Month-year set for Figure C1
sets3 <- construct.sets(mm = matches$mm3, data = subset)

# Sets for DID tables C.1 and C.2
sets <- construct.sets.year(mm = matches$mm3, data_month = subset, data_year = full, years_before = 1, years_after = 3)



#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Appendix: Balance  #### 
# (Table B.3)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

## Treatment/Control sets ----

# construct treatment set
all_treats <- unique(full$sc[full$treat==1])
treats <- unique(full$sc[full$treat==1])
full$year_signed <- as.numeric(substr(full$date1, 1, 4))

t0 <- NULL
for (i in c(1:length(treats))) {
  t <- min(full$year[full$sc==treats[i] & full$treat==1])
  start <- t - 10
  end <- t + 10
  temp <- full[full$sc==treats[i] & full$year>=start & full$year<=end,]
  temp <- temp[order(temp$year),]
  temp$t <- temp$year - t
  t0 <- rbind(t0, temp)
}

# construct control set
controls <- unique(full$sc[full$control==1])
# get rid of two controls that ended up being treated later
controls <- controls[!controls %in% all_treats]

c1 <- NULL
for (i in c(1:length(controls))) {
  t <- min(full$year[full$sc==controls[i] & full$denied==1])
  start <- t - 10
  end <- t + 10
  temp <- full[full$sc==controls[i] & full$year>=start & full$year<=end,]
  temp <- temp[order(temp$year),]
  temp$t <- temp$year - t
  c1 <- rbind(c1, temp)
}

# matched control set
c2 <- sets$control.set

t0$delta5_foreign_pct <- t0$delta5_foreign * 100
c1$delta5_foreign_pct <- c1$delta5_foreign * 100
c2$delta5_foreign_pct <- c2$delta5_foreign * 100

t0$pct_foreign <- (t0$foreign/t0$pop) * 100
c1$pct_foreign <- (c1$foreign/c1$pop) * 100
c2$pct_foreign <- (c2$foreign/c2$pop) * 100

t0$pct_undoc <- (t0$undoc_county/t0$pop) * 100
c1$pct_undoc <- (c1$undoc_county/c1$pop) * 100
c2$pct_undoc <- (c2$undoc_county/c2$pop) * 100

fileConn <- file("balance.tex")
writeLines(c("{", "\n",
             "\\begin{tabular}{ l c c c}", "\n",
             "& Treatment & Control & Control \\\\", "\n",
             "& & (Denied) & (Matched) \\\\", "\n", 
             "\\hline", "\n",
             "\\multicolumn{3}{l}{\\textbf{Population} } \\\\", "\n",
             "\\hspace{2mm} Total & ", 
             prettyNum(round(mean(t0$pop[t0$t==-1], na.rm = TRUE), digits = 0), big.mark = ","), " & $", 
             prettyNum(round(mean(c1$pop[c1$t==-1], na.rm = TRUE), digits = 0), big.mark = ","), "^{", star(t=t0, c=c1, var = "pop"), "}$ ", " & $", 
             prettyNum(round(mean(c2$pop[c2$t==-1], na.rm = TRUE), digits = 0), big.mark = ","), "^{", star(t=t0, c=c2, var = "pop"), "}$ \\\\ ", 
             "\\hspace{2mm} Percent undocumented (county-level estimate) & ",
             prettyNum(round(mean(t0$pct_undoc[t0$t==-1], na.rm = TRUE), digits = 1), big.mark = ","), " & $",
             prettyNum(round(mean(c1$pct_undoc[c1$t==-1], na.rm = TRUE), digits = 1), big.mark = ","), "^{", star(t=t0, c=c1, var = "pct_undoc"), "}$ ", " & $",
             prettyNum(round(mean(c2$pct_undoc[c2$t==-1], na.rm = TRUE), digits = 1), big.mark = ","), "^{", star(t=t0, c=c2, var = "pct_undoc"), "}$ \\\\ ",
             "\\hspace{2mm} Percent foreign-born & ",
             prettyNum(round(mean(t0$pct_foreign[t0$t==-1], na.rm = TRUE), digits = 1), big.mark = ","), " & $",
             prettyNum(round(mean(c1$pct_foreign[c1$t==-1], na.rm = TRUE), digits = 1), big.mark = ","), "^{", star(t=t0, c=c1, var = "pct_foreign"), "}$ ", " & $",
             prettyNum(round(mean(c2$pct_foreign[c2$t==-1], na.rm = TRUE), digits = 1), big.mark = ","), "^{", star(t=t0, c=c2, var = "pct_foreign"), "}$ \\\\ ",
             "\\hspace{2mm} Percentage point change in foreign-born & ",
             prettyNum(round(mean(t0$delta5_foreign_pct[t0$t==-1], na.rm = TRUE), digits = 1), big.mark = ","), " & $",
             prettyNum(round(mean(c1$delta5_foreign_pct[c1$t==-1], na.rm = TRUE), digits = 1), big.mark = ","), "^{", star(t=t0, c=c1, var = "delta5_foreign_pct"), "}$ ", " & $",
             prettyNum(round(mean(c2$delta5_foreign_pct[c2$t==-1], na.rm = TRUE), digits = 1), big.mark = ","), "^{", star(t=t0, c=c2, var = "delta5_foreign_pct"), "}$ \\\\ ",
             "\\hspace{4mm} population, past 5 years \\\\",
             "\\multicolumn{3}{l}{\\textbf{Immigration Enforcement Outcomes} }  \\\\", "\n",
             "\\hspace{2mm} Detainers & ",
             prettyNum(round(mean(t0$ndetainees[t0$t==-1], na.rm = TRUE), digits = 0), big.mark = ","), " & $",
             prettyNum(round(mean(c1$ndetainees[c1$t==-1], na.rm = TRUE), digits = 0), big.mark = ","), "^{", star(t=t0, c=c1, var = "ndetainees"), "}$  ", " & $",
             prettyNum(round(mean(c2$ndetainees[c2$t==-1], na.rm = TRUE), digits = 0), big.mark = ","), "^{", star(t=t0, c=c2, var = "ndetainees"), "}$ \\\\ ",
             "\\hspace{2mm} Detainers with past criminal charge & ",
             prettyNum(round(mean(t0$ncharged[t0$t==-1], na.rm = TRUE), digits = 1), big.mark = ","), " & $",
             prettyNum(round(mean(c1$ncharged[c1$t==-1], na.rm = TRUE), digits = 1), big.mark = ","), "^{", star(t=t0, c=c1, var = "ncharged"), "}$ ", " & $",
             prettyNum(round(mean(c2$ncharged[c2$t==-1], na.rm = TRUE), digits = 1), big.mark = ","), "^{", star(t=t0, c=c2, var = "ncharged"), "}$ \\\\ ",
             "\\multicolumn{3}{l}{\\textbf{Other County Characteristics} }  \\\\", "\n",
             "\\hspace{2mm} Had 2006 immigration protests & ",
             prettyNum(round(mean(t0$prot2006[t0$t==-1], na.rm = TRUE), digits = 2), big.mark = ","), " & $",
             prettyNum(round(mean(c1$prot2006[c1$t==-1], na.rm = TRUE), digits = 2), big.mark = ","), "^{", star(t=t0, c=c1, var = "prot2006"), "}$ ", " & $",
             prettyNum(round(mean(c2$prot2006[c2$t==-1], na.rm = TRUE), digits = 2), big.mark = ","), "^{", star(t=t0, c=c2, var = "prot2006"), "}$ \\\\ ",
             "\\hspace{2mm} Participated in Secure Communities & ",
             prettyNum(round(mean(t0$sec_comm[t0$t==-1], na.rm = TRUE), digits = 2), big.mark = ","), " & $",
             prettyNum(round(mean(c1$sec_comm[c1$t==-1], na.rm = TRUE), digits = 2), big.mark = ","), "^{", star(t=t0, c=c1, var = "sec_comm"), "}$  ", " & $",
             prettyNum(round(mean(c2$sec_comm[c2$t==-1], na.rm = TRUE), digits = 2), big.mark = ","), "^{", star(t=t0, c=c2, var = "sec_comm"), "}$ \\\\ ",
             "\\hspace{2mm} Democratic vote share (\\%), last presidential & ",
             prettyNum(round(mean(t0$pct_dem[t0$t==-1], na.rm = TRUE), digits = 1), big.mark = ","), " & $",
             prettyNum(round(mean(c1$pct_dem[c1$t==-1], na.rm = TRUE), digits = 1), big.mark = ","), "^{", star(t=t0, c=c1, var = "pct_dem"), "}$ ", " & $",
             prettyNum(round(mean(c2$pct_dem[c2$t==-1], na.rm = TRUE), digits = 1), big.mark = ","), "^{", star(t=t0, c=c2, var = "pct_dem"), "}$ \\\\ ",
             "\\hspace{4mm} election & \\\\ ", 
             "\\hline ", "\n", 
             "\\end{tabular}", "\n",
             "}"), 
           sep="", fileConn)
close(fileConn)



#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Appendix Figure B.2: Map of counties in sample  ####
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

county <- full[,!names(full) %in% "sc"]

c1 <- map_data("county")
data(county.regions)

names(county.regions)[c(1,3,4)] <- c("fips", "subregion", "region")
mapdata <- left_join(c1, county.regions, by=c("region", "subregion"))
names(mapdata)[7] <- "GEO_id"

# treated units
treats <- unique(county$GEO_id[county$treat==1])

# controls 
controls <- unique(county$GEO_id[county$denied==1])
controls <- controls[!controls %in% treats]

# entered under trump 
county$year_signed <- as.numeric(substr(county$date1, 1, 4))
trump <- unique(county$GEO_id[county$signer==1 & county$year_signed>=2017])

mapdata$countyg <- 0
mapdata$countyg[mapdata$GEO_id %in% treats] <- 1
mapdata$countyg[mapdata$GEO_id %in% controls] <- 2
mapdata$countyg[mapdata$GEO_id %in% trump] <- 3


## Figure B.2 ----

pdf("map.pdf", width = 8, height = 4)
ggplot(data = mapdata) + 
  geom_polygon(aes(x = long, y = lat, fill = as.factor(countyg), group = group), color = NA, lwd=.1) +
  coord_fixed(1.3) +
  guides(fill=FALSE) +
  theme_void() +
  borders("county", colour = "gray30", size=.2) + 
  borders("state", colour = "gray15", size=.2) + 
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank()) +
  scale_fill_manual(values = c("1" = "red", "2" = "yellow", "3" = "darkorange"))
dev.off()



#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Appendix: Histogram  #### 
# (Figure B.1)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

## plot of raw detainer data
full_my$transferred <- full_my$ice_custody
full_my$not_transferred <- full_my$ndetainees
full_my$date2 <- as.Date(as.yearmon(as.character(full_my$my)))

toplot <- full_my %>% group_by(date2) %>% 
  summarise(transferred = sum(transferred), 
            not_transferred = sum(not_transferred))

brks <- as.Date(c("2002-01-01", "2004-01-01", "2006-01-01", "2008-01-01", "2010-01-01", "2012-01-01", 
                  "2014-01-01", "2016-01-01"))
lbls <- lubridate::year(brks)

pdf("detpermo.pdf", width = 8, height = 5)
ggplot(toplot, aes(x = date2)) + 
  geom_area(aes(y = transferred + not_transferred, fill = "transferred")) + 
  geom_area(aes(y = not_transferred, fill = "not_transferred")) +
  labs(y = "Number of detainer requests",
       x = "Year") + 
  scale_x_date(breaks = brks, labels = lbls) + 
  scale_fill_manual(name = "Transferred to ICE Custody", 
                    labels = c("transferred" = "Yes", "not_transferred" = "No", "na" = "No Data"),
                    values = c("transferred" = "red4", "not_transferred" = "gray40", "na" = "black")) +  
  theme_minimal() + 
  theme(legend.position = "bottom", panel.grid.minor = element_blank())
dev.off()


#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Appendix: Month-year analyses ####
# (Figure C.1, Table C.1, Table C.2)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

## Figure C.1 ----
# Make pretreatment trends figures with matched sets

pdf("pretrends_matched.pdf", width = 5.5, height = 4.2)
trendplot_my(byvar = "log_detainees", 
          treat.set = sets3$treat.set, 
          control.set = sets3$control.set, 
          labels = c("Matched Controls", "Treated Counties"),
          ylab = "Logged detainers")
dev.off()


# Make pretreatment trends figures with matched sets

# overall outcomes
est_all <- compute.did(mm.obj = matches$mm3, byvar = "log_detainees", sets = sets, time.id = "period", unit.id = "id", treat = "treat", data = subset, lag = 1, lead = c(0:3))

# outcomes by facility type 
est_county <- compute.did(mm.obj = matches$mm3, byvar = "log_detainees_county", sets = sets, time.id = "period", unit.id = "id", treat = "treat", data = full, lag = 1, lead = c(0:3))
est_noncounty <- compute.did(mm.obj = matches$mm3, byvar = "log_detainees_noncounty", sets = sets, time.id = "period", unit.id = "id", treat = "treat", data = full, lag = 1, lead = c(0:3))

# past criminal record
est_nocharge <- compute.did(mm.obj = matches$mm3, byvar = "log_notcharged_county", sets = sets, time.id = "period", unit.id = "id", treat = "treat", data = full, lag = 1, lead = c(0:3))
est_ncharged_low <- compute.did(mm.obj = matches$mm3, byvar = "log_ncharged_low_county", sets = sets, time.id = "period", unit.id = "id", treat = "treat", data = full, lag = 1, lead = c(0:3))
est_ncharged_med <- compute.did(mm.obj = matches$mm3, byvar = "log_ncharged_med_county", sets = sets, time.id = "period", unit.id = "id", treat = "treat", data = full, lag = 1, lead = c(0:3))
est_ncharged_high <- compute.did(mm.obj = matches$mm3, byvar = "log_ncharged_high_county", sets = sets, time.id = "period", unit.id = "id", treat = "treat", data = full, lag = 1, lead = c(0:3))

# Table C.1
fileConn <- file("did_main.tex")
writeLines(c("{", "\n",
             "\\begin{tabular}{ l c c c c }", "\n",
             "& $t$ & $t+1$ & $t+2$ & $t+3$ \\\\", "\n", 
             "\\hline", "\n", 
             "\\\\", "\n", 
             "Logged detainers & ", 
             "$", format(round(est_all$pe["t+0", "estimate"], 1), nsmall = 2), "^{", est_all$pe["t+0", "star"], "}$ & ", 
             "$", format(round(est_all$pe["t+1", "estimate"], 1), nsmall = 2), "^{", est_all$pe["t+1", "star"], "}$ & ", 
             "$", format(round(est_all$pe["t+2", "estimate"], 1), nsmall = 2), "^{", est_all$pe["t+2", "star"], "}$ & ", 
             "$", format(round(est_all$pe["t+3", "estimate"], 1), nsmall = 2), "^{", est_all$pe["t+3", "star"], "}$ ", " \\\\ ", "\n", 
             "& \\small (", format(round(est_all$pe["t+0", "std.error"], 1), nsmall = 2), 
             ") & \\small (", format(round(est_all$pe["t+1", "std.error"], 1), nsmall = 2), 
             ") & \\small (", format(round(est_all$pe["t+2", "std.error"], 1), nsmall = 2), 
             ") & \\small (", format(round(est_all$pe["t+3", "std.error"], 1), nsmall = 2), 
             ") \\\\", "\n",  					   
             "\\\\", 
             "Logged detainers, county facilities & ", 
             "$", format(round(est_county$pe["t+0", "estimate"], 1), nsmall = 2), "^{", est_county$pe["t+0", "star"], "}$ & ", 
             "$", format(round(est_county$pe["t+1", "estimate"], 1), nsmall = 2), "^{", est_county$pe["t+1", "star"], "}$ & ", 
             "$", format(round(est_county$pe["t+2", "estimate"], 1), nsmall = 2), "^{", est_county$pe["t+2", "star"], "}$ & ", 
             "$", format(round(est_county$pe["t+3", "estimate"], 1), nsmall = 2), "^{", est_county$pe["t+3", "star"], "}$ ", " \\\\ ", "\n", 
             "& \\small (", format(round(est_county$pe["t+0", "std.error"], 1), nsmall = 2), 
             ") & \\small (", format(round(est_county$pe["t+1", "std.error"], 1), nsmall = 2), 
             ") & \\small (", format(round(est_county$pe["t+2", "std.error"], 1), nsmall = 2), 
             ") & \\small (", format(round(est_county$pe["t+3", "std.error"], 1), nsmall = 2), 
             ") \\\\", "\n",  					   
             "\\\\", 
             "Logged detainers, non-county facilities & ", 
             "$", format(round(est_noncounty$pe["t+0", "estimate"], 1), nsmall = 2), "^{", est_noncounty$pe["t+0", "star"], "}$ & ", 
             "$", format(round(est_noncounty$pe["t+1", "estimate"], 1), nsmall = 2), "^{", est_noncounty$pe["t+1", "star"], "}$ & ", 
             "$", format(round(est_noncounty$pe["t+2", "estimate"], 1), nsmall = 2), "^{", est_noncounty$pe["t+2", "star"], "}$ & ", 
             "$", format(round(est_noncounty$pe["t+3", "estimate"], 1), nsmall = 2), "^{", est_noncounty$pe["t+3", "star"], "}$ ", " \\\\ ", "\n", 
             "& \\small (", format(round(est_noncounty$pe["t+0", "std.error"], 1), nsmall = 2), 
             ") & \\small (", format(round(est_noncounty$pe["t+1", "std.error"], 1), nsmall = 2), 
             ") & \\small (", format(round(est_noncounty$pe["t+2", "std.error"], 1), nsmall = 2), 
             ") & \\small (", format(round(est_noncounty$pe["t+3", "std.error"], 1), nsmall = 2), 
             ") \\\\", "\n",  					   
             "\\\\", 
             "\\end{tabular}", "\n",
             "}"), 
           sep="", fileConn)
close(fileConn)

# Table C.2
fileConn <- file("did_charge.tex")
writeLines(c("{", "\n",
             "\\begin{tabular}{ l c c c c }", "\n",
             "& $t$ & $t+1$ & $t+2$ & $t+3$ \\\\", "\n", 
             "\\hline", "\n", 
             "\\\\", "\n", 
             "Logged detainers, no past criminal charge & ", 
             "$", format(round(est_nocharge$pe["t+0", "estimate"], 2), nsmall = 2), "^{", est_nocharge$pe["t+0", "star"], "}$ & ", 
             "$", format(round(est_nocharge$pe["t+1", "estimate"], 2), nsmall = 2), "^{", est_nocharge$pe["t+1", "star"], "}$ & ", 
             "$", format(round(est_nocharge$pe["t+2", "estimate"], 2), nsmall = 2), "^{", est_nocharge$pe["t+2", "star"], "}$ & ", 
             "$", format(round(est_nocharge$pe["t+3", "estimate"], 2), nsmall = 2), "^{", est_nocharge$pe["t+3", "star"], "}$ ", " \\\\ ", "\n", 
             "& \\small (", format(round(est_nocharge$pe["t+0", "std.error"], 2), nsmall = 2), 
             ") & \\small (", format(round(est_nocharge$pe["t+1", "std.error"], 2), nsmall = 2), 
             ") & \\small (", format(round(est_nocharge$pe["t+2", "std.error"], 2), nsmall = 2), 
             ") & \\small (", format(round(est_nocharge$pe["t+3", "std.error"], 2), nsmall = 2), 
             ") \\\\", "\n",  					   
             "\\\\", 
             "Logged detainers, past criminal charge, misdemeanor & ", 
             "$", format(round(est_ncharged_low$pe["t+0", "estimate"], 2), nsmall = 2), "^{", est_ncharged_low$pe["t+0", "star"], "}$ & ", 
             "$", format(round(est_ncharged_low$pe["t+1", "estimate"], 2), nsmall = 2), "^{", est_ncharged_low$pe["t+1", "star"], "}$ & ", 
             "$", format(round(est_ncharged_low$pe["t+2", "estimate"], 2), nsmall = 2), "^{", est_ncharged_low$pe["t+2", "star"], "}$ & ", 
             "$", format(round(est_ncharged_low$pe["t+3", "estimate"], 2), nsmall = 2), "^{", est_ncharged_low$pe["t+3", "star"], "}$ ", " \\\\ ", "\n", 
             "& \\small (", format(round(est_ncharged_low$pe["t+0", "std.error"], 2), nsmall = 2), 
             ") & \\small (", format(round(est_ncharged_low$pe["t+1", "std.error"], 2), nsmall = 2), 
             ") & \\small (", format(round(est_ncharged_low$pe["t+2", "std.error"], 2), nsmall = 2), 
             ") & \\small (", format(round(est_ncharged_low$pe["t+3", "std.error"], 2), nsmall = 2), 
             ") \\\\", "\n",  					   
             "\\\\", 
             "Logged detainers, past criminal charge, felony (low) & ", 
             "$", format(round(est_ncharged_med$pe["t+0", "estimate"], 2), nsmall = 2), "^{", est_ncharged_med$pe["t+0", "star"], "}$ & ", 
             "$", format(round(est_ncharged_med$pe["t+1", "estimate"], 2), nsmall = 2), "^{", est_ncharged_med$pe["t+1", "star"], "}$ & ", 
             "$", format(round(est_ncharged_med$pe["t+2", "estimate"], 2), nsmall = 2), "^{", est_ncharged_med$pe["t+2", "star"], "}$ & ", 
             "$", format(round(est_ncharged_med$pe["t+3", "estimate"], 2), nsmall = 2), "^{", est_ncharged_med$pe["t+3", "star"], "}$ ", " \\\\ ", "\n", 
             "& \\small (", format(round(est_ncharged_med$pe["t+0", "std.error"], 2), nsmall = 2), 
             ") & \\small (", format(round(est_ncharged_med$pe["t+1", "std.error"], 2), nsmall = 2), 
             ") & \\small (", format(round(est_ncharged_med$pe["t+2", "std.error"], 2), nsmall = 2), 
             ") & \\small (", format(round(est_ncharged_med$pe["t+3", "std.error"], 2), nsmall = 2), 
             ") \\\\", "\n",  					   
             "\\\\", 
             "Logged detainers, past criminal charge, felony (high) & ", 
             "$", format(round(est_ncharged_high$pe["t+0", "estimate"], 2), nsmall = 2), "^{", est_ncharged_high$pe["t+0", "star"], "}$ & ", 
             "$", format(round(est_ncharged_high$pe["t+1", "estimate"], 2), nsmall = 2), "^{", est_ncharged_high$pe["t+1", "star"], "}$ & ", 
             "$", format(round(est_ncharged_high$pe["t+2", "estimate"], 2), nsmall = 2), "^{", est_ncharged_high$pe["t+2", "star"], "}$ & ", 
             "$", format(round(est_ncharged_high$pe["t+3", "estimate"], 2), nsmall = 2), "^{", est_ncharged_high$pe["t+3", "star"], "}$ ", " \\\\ ", "\n", 
             "& \\small (", format(round(est_ncharged_high$pe["t+0", "std.error"], 2), nsmall = 2), 
             ") & \\small (", format(round(est_ncharged_high$pe["t+1", "std.error"], 2), nsmall = 2), 
             ") & \\small (", format(round(est_ncharged_high$pe["t+2", "std.error"], 2), nsmall = 2), 
             ") & \\small (", format(round(est_ncharged_high$pe["t+3", "std.error"], 2), nsmall = 2), 
             ") \\\\", "\n",  					   
             "\\\\", 
             "\\end{tabular}", "\n",
             "}"), 
           sep="", fileConn)
close(fileConn)


#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Appendix: Print treated and control units
# For Table B.1 and B.2
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# Table B.1
unique(sets3$treat.set[,c("sc", "date1")])

# Table B.2
b2_temp <- unique(full[full$control==1 & !full$GEO_id %in% treats & full$denied==1, c("sc", "year")])
b2 <- b2_temp %>%
  group_by(sc) %>%
  filter(year == min(year, na.rm = TRUE)) %>%
  ungroup()

